#include "solver.h"

#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <stdio.h>

#include "portable_time.h"

double* solve_withoutOMP(double* u0, Area area, double epsilon, TwoDimFunc f)
{
    int n = area.n;
    double a = area.a, b = area.b;
    double* u_cur = new double[n * n];
    double* u_prev = new double[n * n];
    double h = (b - a) / n;
    double error;
    double start, finish;

    double* f_cur = new double[n * n];

    start = PortableGetTime();
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            f_cur[i * n + j] = f(a + i * h, a + j * h);
        }
    }
    finish = PortableGetTime();
    printf("Compute function time: %.4f\n", finish - start);

    memcpy(u_prev, u0, n * n * sizeof(double));
    memcpy(u_cur, u0, n * n * sizeof(double));

    start = PortableGetTime();
    do
    { 
        for (int i = 1; i < n - 1; ++i)
        {
            for (int j = 1; j < n - 1; ++j)
            {
                u_cur[i * n + j] = 0.25 * 
                    (u_prev[(i - 1) * n + j] + u_prev[(i + 1) * n + j] +
                    u_prev[i * n + j - 1] + u_prev[i * n + j + 1] - 
                    h * h * f_cur[i * n + j]);
            }
        }

        double e;
        error = abs(u_cur[0] - u_prev[0]);
        for (int i = 1; i < n * n - 1; ++i)
        {
            e = abs(u_cur[i] - u_prev[i]);
            if (e > error) error = e;
        }

        double* tmp;
        tmp = u_cur; u_cur = u_prev; u_prev = tmp;

    } while (error > epsilon);
    finish = PortableGetTime();
    printf("Solve time: %.4f\n", finish - start);

    delete[] u_cur;
    return u_prev;
}